library(topGO)
library(dplyr)
library(tidyr)
library(RColorBrewer)
library(ClassDiscovery)
# read in all the files
setwd("~/Desktop/brainRNASeq/R_analysis/R_output/topGO/SCN")
BP.GO.data <- readRDS("BP/SCN_all_comparsion_topGO_BP.RData")
KS.result <- read.csv("BP/SCN_all_comparsion_BP_KS.csv", header = TRUE, stringsAsFactors = FALSE)
# read in the rlog transformation table
setwd("~/Desktop/brainRNASeq/R_analysis/R_output/results_tables/")
rlog.transformed.data <- read.csv("SCN/SCN_rlog_transformated_data.csv", header = TRUE, 
    stringsAsFactors = FALSE)
# print the result table
KS.result
##         GO.ID                                  Term Annotated Significant
## 1  GO:0051641                 cellular localization      1869         384
## 2  GO:0046907               intracellular transport      1135         229
## 3  GO:0033036            macromolecule localization      1997         377
## 4  GO:0051649 establishment of localization in cell      1402         294
## 5  GO:0008104                  protein localization      1811         347
## 6  GO:0045184 establishment of protein localization      1283         247
## 7  GO:0031175         neuron projection development       804         167
## 8  GO:0071702           organic substance transport      1631         298
## 9  GO:0048666                    neuron development       888         179
## 10 GO:0071705           nitrogen compound transport      1426         268
##    Expected Rank.in.classicKS classicKS classicFisher elimFisher elimKS
## 1    356.34                 1   7.7e-11         0.040       0.75 0.1660
## 2    216.40                 2   8.5e-11         0.167       0.24 0.0024
## 3    380.75                 3   8.3e-10         0.604       0.93 0.1761
## 4    267.30                 4   1.3e-09         0.029       0.28 0.8363
## 5    345.28                 5   1.3e-09         0.466       0.89 0.3271
## 6    244.62                 6   4.3e-09         0.441       0.64 0.0436
## 7    153.29                 7   9.0e-09         0.109       0.15 0.0105
## 8    310.97                 8   9.6e-09         0.822       0.92 0.0231
## 9    169.31                 9   9.9e-09         0.205       0.27 0.0075
## 10   271.88                10   1.1e-08         0.623       0.80 0.2805
# extract the list of genes by the topGO

gene.from.GOterms <- genesInTerm(BP.GO.data, KS.result$GO.ID)

colnames(rlog.transformed.data)[1] <- "ensembl_gene_id"
for (i in 1:10) {
    # looping for each GO accessing the list of gene in each GO term
    gene.in.GO <- as.data.frame(gene.from.GOterms[i])
    gene.in.GO <- as.vector(gene.in.GO[, 1])
    
    # subset the genes
    gene.interested.rlog <- rlog.transformed.data[rlog.transformed.data$ensembl_gene_id %in% 
        gene.in.GO, ]
    
    if (length(gene.interested.rlog$ensembl_gene_id) >= 2) {
        # convert the dataframe to matrix with only the genes in the GO term
        rownames(gene.interested.rlog) <- gene.interested.rlog$ensembl_gene_id
        gene.interested.rlog <- select(gene.interested.rlog, -ensembl_gene_id)
        gene.interested.rlog <- as.matrix(gene.interested.rlog)
        
        # ploting the heatmap
        aspectHeatmap(gene.interested.rlog, scale = "row", Colv = NA, hExp = 1, 
            wExp = 1.5, margins = c(6, 10), main = paste(KS.result$GO.ID[i]), 
            col = colorRampPalette(rev(brewer.pal(9, "RdBu")))(255))
    }
}

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ClassDiscovery_3.3.12 oompaBase_3.2.9       cluster_2.1.0        
##  [4] RColorBrewer_1.1-2    tidyr_1.0.0           dplyr_0.8.3          
##  [7] topGO_2.36.0          SparseM_1.77          GO.db_3.8.2          
## [10] AnnotationDbi_1.46.1  IRanges_2.18.2        S4Vectors_0.22.1     
## [13] Biobase_2.44.0        graph_1.62.0          BiocGenerics_0.30.0  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2         oompaData_3.1.1    compiler_3.6.0    
##  [4] pillar_1.4.2       formatR_1.7        tools_3.6.0       
##  [7] mclust_5.4.5       zeallot_0.1.0      digest_0.6.21     
## [10] bit_1.1-14         lifecycle_0.1.0    RSQLite_2.1.2     
## [13] evaluate_0.14      memoise_1.1.0      tibble_2.1.3      
## [16] lattice_0.20-38    pkgconfig_2.0.3    rlang_0.4.0       
## [19] DBI_1.0.0          yaml_2.2.0         xfun_0.9          
## [22] stringr_1.4.0      knitr_1.25         vctrs_0.2.0       
## [25] tidyselect_0.2.5   bit64_0.9-7        grid_3.6.0        
## [28] glue_1.3.1         R6_2.4.0           rmarkdown_1.15    
## [31] purrr_0.3.2        blob_1.2.0         magrittr_1.5      
## [34] matrixStats_0.55.0 backports_1.1.4    htmltools_0.3.6   
## [37] assertthat_0.2.1   stringi_1.4.3      crayon_1.3.4